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WHEN INTERPOLATION-INDUCED REFLECTION ARTIFACT MEETS 
TIME-FREQUENCY ANALYSIS 

YU-TING LIN, PATRICK FLANDRIN, AND HAU-TIENG WU 


Ab STRACT. While extracting the temporal dynamical features based on the time-frequency 
analyses, like the reassignment and synchrosqueezing transform, attracts more and more 
interest in bio-medical data analysis, we should be careful about artifacts generated by in¬ 
terpolation schemes, in particular when the sampling rate is not significantly higher than 
the frequency of the oscillatory component we are interested in. In this study, we for¬ 
mulate the problem called the reflection effect and provide a theoretical justification of 
the statement. We also show examples in the anesthetic depth analysis with clear but un¬ 
desirable artifacts. The results show that the artifact associated with the reflection effect 
exists not only theoretically but practically. Its influence is pronounced when we apply the 
time-frequency analyses to extract the time-varying dynamics hidden inside the signal. In 
conclusion, we have to carefully deal with the artifact associated with the reflection effect 
by choosing a proper interpolation scheme. 


1. Introduction 

It has been widely accepted that several aspects of the health status could be well ob¬ 
served by analyzing recorded physiological time series. In particular, the time-varying 
oscillatory pattern inside the electrocardiogram (ECG) or respiratory signal contains abun¬ 
dant health information, for example, the heart rate variability (HRV) El El 13 hidden in¬ 
side the R peak to R peak interval (RRI) time series and the instantaneous heart rate (IHR), 
the breathing pattern variability (BPV) representing the time varying rate of the respira¬ 
tory signal El Elia. It is well known that power spectrum is not a suitable tool when 
the time-varying dynamics in the signal is the main target to analyze, as power spectrum 
reflects only the global oscillatory information, and hence could not properly extract the 
dynamical information, which is local in nature. In general, a popular and powerful way 
to study the time-varying oscillatory pattern inside a time series is the time-frequency (TF) 
analysis, which allows us to efficiently extract how a signal oscillates at each time instant. 
There have been several TF analysis techniques proposed, including linear methods like 
short time Fourier transform (STFT), continuous wavelet transform (CWT) □13 and mul¬ 
tiwindow approach a, the quadratic methods like Wigner-Ville distribution and Cohen 
class m, nonlinear methods like reassignment (RM) technique oni[n]i, synchrosqueezing 
transform (SST) [TJKTSl, multi-tapered RM ifllll . multi-tapered SST |[T3[T6l, ConceFT 
EtI . empirical model decomposition E3, sparse TF analysis CD, iterative filtering |[20l . 
etc. The potential of the TF analysis has been shown in several different fields, in particular 
the medical field; for example, EH Ea 123 da, to name but a few. 

While there exists a lot of information in the physiological signals we could easily 
approach, in some signals, like IHR, there are two particular features that should not be 
neglected. First, they are sampled in a non-uniform fashion; second, in many situations, 
they are often sampled at a rate which is not significantly high. Typical examples include 
the IHR and the ECG-derived respiratory (EDR) signal extracted from the ECG signal, the 
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IHR estimated from the photoplethysmography signal 1241 . etc, in which case the sam¬ 
pling is non-uniform and the sampling rate is determined by the heart rate. In order to 
apply the TF analysis to study these observed signals, a common practice is to apply the 
digital-to-analogue conversion to recover the original continuous signal EH [2 EH EH na, 
for example, the spline interpolation. Among different choices of the spline interpolation 
schemes, the cubic spline interpolation is commonly chosen as it balances between the in¬ 
terpolation property and the overfitting issue raised by the interpolation. On the other hand, 
the interpolation is commonly associated with a convolutional kernel. While it is generally 
accepted that the quality of recovering the underlying signal is good by a suitably chosen 
interpolation scheme, it might cause some undesirable artifacts, for example, by the side 
lobe effect inherited in the spline interpolation scheme |E21- Thus, although the TF analy¬ 
sis have been successfully applied to several physiological problems for different purposes, 
we should be careful about the analysis results. For example, while the main purpose of 
the nonlinear TF analysis techniques like RM and SST is to sharpen/enhance the TF rep¬ 
resentation determined by the linear TF analysis, these techniques might also enhance the 
artifact generated by the interpolation scheme and hence mislead the interpretation. This 
problem is further complicated by the possible non-uniform sampling scheme. 

To demonstrate the potential problem caused by this issue, here we show a real example 
in the respiratory signal analysis. The airfiow respiratory signal and the ECG signal are 
recorded simultaneously, and we obtain the FDR signal from the ECG signal. The FDR 
signal is generated by interpolating the amplitudes of detected R peaks via the cubic spline 
interpolation, and is influenced by the inevitable sampling effect inherited in the R peak 
location, which is not only non-uniform but also with low sampling rate. In Figure the 
respiratory signal and the FDR signal are shown together, as well as their TF representa¬ 
tions determined by the multi-tapered RM. Clearly, while these two signals “look” similar 
in the sense of the “fast-slow” oscillatory pattern, their TF representations are very differ¬ 
ent. While the TF representation of the respiratory signals show the multiples of the base 
respiratory frequency at about 0.5Hz, the TF representation of the FDR signal show “two 
different components”, where the component with the higher frequency is not the multiple 
of the base respiratory frequency at about 0.5Hz. Clearly, the TF representation of the FDR 
signal will mislead the interpretation. 

Thus, to extract the correct information from this kind of time series, it is essential 
to understand the infiuence on the analysis results caused by the sampling scheme, the 
interpolation scheme and the TF analysis. In this paper, the refiection effect caused by the 
interpolation scheme is formalized, how the sampling scheme is involved is discussed, and 
a theoretical justification of this effect is provided. In addition, we show several medical 
examples where the artifacts might mask the whole interpretation. 

This paper is organized in the following way. In Sectionj^ the adaptive harmonic model 
is introduced to model the commonly encountered oscillatory signals. The notion instan¬ 
taneous Nyquist rate is introduced to quantify the nonuniform sampling that serves as a 
framework for our analysis. In Sectionthe refiection effect is formalized and a theoret¬ 
ical justification is provided. In addition, a solution to this artificial effect is proposed. A 
series of numerical evidences, as well as real examples from the anesthesia are provided in 
Sections]^ and [H In SectionE] we conclude the paper with a series of discussions. 

2. The adaptive harmonic model and instantaneous Nyquist rate 

Among different features of time series, oscillatory pattern is the main target of several 
TF analysis techniques. When it comes to the oscillation, Fourier analysis is commonly 
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Figure 1. The time-frequency representations of the ECG derived res¬ 
piratory (EDR) signal and the airflow signal determined by the multi¬ 
taper reassignment (RM). The signals are recorded simultaneously. Top: 
the EDR signal based on the cubic spline interpolation from the R peaks 
amplitudes determined from the lead II ECG signal. The mean of the 
EDR signal is removed. Middle top: the multi-taper RM of the EDR. 

The instantaneous Nyquist frequency (INE) is superimposed as a red 
dashed curve. Middle bottom: the airflow signal recorded simultane¬ 
ously. At around 150 second, an automatic calibration happens, which 
leads to a short zero period. Bottom: the multi-taper RM of the airflow 
signal. Note that the breathing rate is below 0.5 Hz in the beginning and 
increases gradually. It is clear that the base component with instanta¬ 
neous frequency (IE) about 0.5 Hz in the airflow signal is well captured 
by the EDR signal, while in the EDR signal there is an artiflcial compo¬ 
nent with IE higher than INE, which mirrors the base component IE via 
the INE. 

the first choice among others. However, in general, the oscillatory pattern might change 

according to time, and its momentary behavior might not be easily captured by the Eourier 
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technique and new techniques are needed. While the time-varying oscillatory pattern might 
be as complicated as one could imagine, in medical applications it is commonly controlled 
by the physiological constraints. For example, the heart rate could not be as fast as possible, 
and the amplitude of the respiratory signal is limited by the lung capacity. Thus, to quantify 
the oscillatory pattern inside an observed time series, we consider the following adaptive 
harmonic model ifT^ [TSl : 

Definition 2.1. Fix 0 < e <C 1, £ ^ ci < C 2 . The functional class, c (M) nL°°(M), 
contains functions of the form a(r)cos(2;r0(r)), where a G C^{R) r\L°°{R), (j) G C^(M), 
Cl < a(t) < C 2 , Cl < 0'(^) < C 2 W{t) \ < and \(l)'\t) \ < for all time ^ G M. 

For /(r) = a{t) cos{27r(j){t)) G we call the instantaneous frequency (IF) and 

a{t) the amplitude modulation (AM) of f{t). In the adaptive harmonic model, locally all 
functions in behave like a harmonic function (or sine wave), and the deviation from 

being a harmonic function is controlled by £. We call a function in an intrinsic mode 

type (IMT) function. A more general model consisting of functions with multiple IMT 
functions and more theoretical discussions could be found in CSlini; the identifiability 
issue of the adaptive harmonic model and the noise and trend model issue have been studied 
and discussed in l(T^ . While the following discussion carries over into to the more general 
model, in this paper, to simplify the discussion, we focus on , where there is only 

one IMF inside the signal. The adaptive harmonic model has been applied to study the 
“non-stationary” physiological dynamics lElElElllaSEllISl. Before proceeding, we 
comment that most time series we could acquire in the real world are real, so we consider 
only real model but not the exponential complex model. We mention that obtaining the 
imaginary part of a given real signal to recover its exponential complex form is not a trivial 
problem ISOllSTlI^ . 

Next we discuss the sampling effect. Take f{t) = a(t) cos{271 (j){t)) G Consider a 

monotonic increasing function y/ eC^ (M). Define a sequence of sampling points {tm}mez 
so thaX tm = il/~^ (m). With this sampling scheme, we obtain samples ^ := {tm,f{im)}mez- 
Note that if \i/{t) = kt, where k > 0, then the sampling is uniform and we get a sample every 
1/k second. To study how much information about / we could obtain from the sampled 
dataset, when i/r is a linear function and / is a band-limited function, we could consider 
the Nyquist-Shannon theory. However, the application of Nyquist-Shannon theory is not 
efficient in our case, since a generic function in has a non-compact supporj^in the 

Fourier domain. The application of the Nyquist-Shannon theory is more difficult in the 
case when y/ is nonlinear; that is, when the sampling is non-uniform. Thus, we consider 
another definition to describe the sampling scheme which reflects the momentary nature of 
dynamical analysis. 

Definition 2.2. Fix c > 0 and 0 < e ^ c. Take y/ G C^{R) so that c < y/'{t) and | y/"{t) \ < 
£y/'{t) for all ^ G M. We call y/'{t) the instantaneous sampling rate (ISR) and (t)j2 the 
instantaneous Nyquist frequency (INF). 

Note that ISR and INF naturally generalize the notion of sampling rate and Nyquist 
frequency - the higher the ISR is at a moment, the higher the sampling rate is around 
this moment. The condition \yf"{t)\ <£yf'{t) says that locally the sampling is close to a 
uniform sampling. A natural problem regarding the above definition of ISR and INF is 
the identifiability issue. Precisely, given a sampling points {tm}m^z^ we could find many 
different functions y/ which leads to the same sampling points. We claim that under the 
provided condition of ISR and INF, they are well-defined up to an error of order £. 

%or the definition of support, we refer the reader to EH P. 284], 
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Theorem 2.1. Fix c > 0. If i//, i// G C^(]R) both satisfy the conditions for ISR and generate 
the same sampling points then we have | y/'it) — ^'{t) \ < 2e and | — \j/{t) \ < 

2£{tm+\ — tm) for all time t G ^m+i] • Globally, we have \ \if'{t) — {t) \ < 2e and \ \i/{t) — 

<2£lc for all time t eR. 


Proof. Note that we have = W{^m) =t for all m G Z. Consider t G By a 

direct calculation, we have 

(1) \w'{t)-^{t)\< f \il/'{s)-^'{s)\ds 

tm 

<e f {x/(s) + {s))ds 

^ tm 

= e[(v^(0 - v{tm)) + {wit) - w{tm))] 

< 2£{w{tm+\) - W{tm)) = 2e, 


where the second inequality holds by the assumption that \ \ < £y/'{t) and \ \ < 

both hold, and the last inequality holds by the monotonic assumption of if/ and 
\j/. To finish the proof, note that by the mean value theorem, we have — y/fm) = 

— tm) for some F G which leads to 


( 2 ) 


tm-\-l tjfi 


¥(tm+l) - ¥(tm) 

¥'(t') 


1 1 
y/'f') ~ c’ 


where the last inequality holds by the assumption that y/'(t) > c and y/'(t) > c. Thus, we 
have 


(3) \w{t)-w{t)\< [ \w’{s)-'‘r{s)\ds 

Jtm 

^ 2£{t tm) ^ 2£{tm-\-l ^m) ^ 2^/ C 


and hence the proof is done. 


□ 


This theorem essentially says that the ISR for a given sampling points {tm}mez is well 
defined up to the error of order £. We mention that for a randomly given sampling points 
{tm}mez SO that tfn^i > Hiay not be able to define a ISR function which satisfies the 

In this paper, we focus on sampling scheme which satisfies 


condition in Definition 2.2 
Definition 12.21 


In addition to ISR and INF, note that we could also naively generalize the notion of 
Nyquist rate to its instantaneous version. 


Definition 2.3. Take a signal f{t) = a{t) cos(2;r0(r)) G We call 2(1)'(t) the instan¬ 

taneous Nyquist rate (INR) of the function f{t). 


Note that the INR reduces to the notion of Nyquist rate when 0(r) is linear and a{t) 
is constant. We will always assume that yf\t) > 2(^'{t)\ that is, at each time instant, we 
have at least two sampling points from an oscillation, otherwise the oscillatory information 
might be lost. Thus, this is again a natural generalization of the sampling theory under the 
uniform sampling scheme and the band-limited assumption. 


3. The reflection effect 

With the adaptive harmonic model and the samples, we would like to study the underly¬ 
ing dynamical features, like the IF and AM of the signal. A common practice to convert the 
discretized sampling to a continuous function is via an interpolation scheme. We now 
show the structured artifacts caused by the commonly applied spline interpolation scheme. 





6 


YU-TING LIN, PATRICK FLANDRIN, AND HAU-TIENG WU 


In particular, the IF of an artificial component in the interpolated signal is a refiection of 
0' associated with the INF. Further, while the INF is closer to the IF, this artifact is more 
severe, which is the case we encounter when we study the IHR or FDR signals. 

Theorem 3.1. Take/(r) =a{t) cos{27r G and sample/(r) with the ISR > 

2(1)'(t). Fix a ^-th order spline interpolation, where n>l, and denote the interpolated signal 
as Then, we have 

(4) /„(0 = a{t) £ Cn{k)cos[2nikxi/it) - <?>(/))] + 0(Ve), 

where Cn{k) G M depends on n. 


Proof. The proof is divided into three steps. First, we show the proof when the signal is 
harmonic and the sampling is uniform; second, we show the result when the signal is in the 
adaptive harmonic model and the sampling is uniform; third, we show the general proof. 

We start from the harmonic signal f{t) = cos(2;rar), where 0 < a < 1/2. Without loss 
of generality, we assume that the sampling rate is 1 Hz; that is, y/f) = t and the ISR is 1. 
If not, for example, y/f) = kt, where k> 1, we could up wrap the time axis hy s = kt, and 
get y/{s) = s and f{s) = cos{27l{a/k)s), and the argument holds. Note that / is a band- 
limited harmonic function which is also in Note that the IF of f{t) is a, which is 

less than y/'{t) 12. The uniform sampling scheme corresponds to the weighted Dirac train, 
which is a tempered distribution f§ := L/ez/(0^/’ where 5/ is the Dirac delta measure 
supported at / G Z. The ^-th order spline interpolation scheme, where ^ G N, is performed 
as a convolution of f§ with the ^-th order fundamental cardinal spline function which 

satisfies 


(5) 


^(0 = 



sin(;r((^ —/)) 


V 


n+\ 


where f/^ means the Fourier transform of rji^y^y See O (2.14)] or |[35l (4.6.9)] for exam¬ 
ple. Precisely, the interpolated signal based on the ^-th order spline interpolation, denoted 
as /, satisfies / = i^f§, where ^ is the convolution. Note that it is an interpolation and 

we have /(/) = /(/) for all / G Z. By a direct calculation, we have fs{^) = Hhez /(^+^) = 

2 LkeZ [4+ct + 4-a], which leads to /(^) = i Zkez ^)(4) [4+a + 4-ct]• As a result, we 
have 


( 6 ) f(t) = a) cos{2n{k - a)t) 

k^'Zi 

We thus finish the proof when the signal is harmonic and the sampling is uniform. Indeed, 
the main reflection component associated with INF, 1 /2, is T7^(l — a) cos(2;r(l — a)t). 
Second, without loss of generality, we consider a general function f{t) = a{t) cos(2;r0 (t)) G 
0^(0 < 1 /^- Again, we consider a uniform sampling scheme at the sampling 
rate 1 Hz; that is, y/{t) = t. Take /q G Z. Consider a local harmonic approximation of / 
around /q, denoted as f^^^\t) = (2(/o)cos[2;r(0(/o) — adaptive 

harmonic model and the same argument as that in ifT^ . we have the control between / and 
f{k) for all 5’ G M 

(7) f{h + s)=P^\h + s)+C\s\e, 

where C depends on C 2 in the definition of and is uniformly bounded for all /q. Also 

recall the fact that the ^-th order cardinal spline decays exponentially (351 (4.6.2)]. By 
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taking the facts of exponential decay of T 7 („) and the Taylor expansion ([^, we know that 
for t G [/q — 1 /2, Zq + 1 /2], 


( 8 ) 


/eZ 


< 2Ce £ ke-‘=’^ = C'e 


for some constant c, C' > 0. Hence, by ^ we have for all t G [/q — 1/2, /q + 1/2) 

fit) = = (T]w */?^)(0 + O(e) 

(9) = Y,tl{f)ik-<l>'ik))a{l)cos[27t{kl-^{t))] + 0{e), 


where we use the fact that a{t{)) cos[ 2 ;r (0 (/q) — /o 0 '(/o) + (^+ 0 '(/o))O] = ^(0 cos[ 2 ;r (0 {t) + 
kt)]^ 0{£). As a result, we have the proof when / G and the sampling is uniform. 

Indeed, the main reflection component associated with INF, 1/2, isf/^(l — 0 '(/o))f 2 (r)cos[ 2 ;r(r — 
0(r))], which has the IF 1 — 

Third, we consider the case when the sampling is non-uniform and the signal satisfles 
the adaptive harmonic model. We consider f{t) = a(t)co^{27t^{t)) G and the ISR 

y/ so that y/'it) > 2(^'(t) for all t. Denote the sampling points as ^ := where 

yf{ti) = i. Without loss of generality, we focus on T so that 0 '(t) < 1/2 and V^'(t) = 1 > 

20'(t); the other cases could be proved by scaling. To simplify the notation, we assume 
that T = A/2, where N = 2[l/v^]. By the assumption of i/r, we know that over the 
interval [0, A], | yf'{s) — 11 < for all 5’ G [0,A]. Denote / = {/ G Z| 0 < < A} C Z and 

/={0,1,...,A} C Z. Note that \ti — l\ <C'\1 — t\£ for all / G /, where C' is a constant 
depending on i/r, so j^/ — / | is bounded by y/e and we know that |/| = |/| = A+1. Denote 
Nnj to be the B-spline of order n deflned on ^ (361 (10.2.32)]; that is. 


(r — tiY 

(1^) '•= -tj) ^ j+^+1 . ^ X ’ 

k=j ^k^l=jVl-^k) 

where x\ is the truncated power deflned as x\ = x^ when a > 0 and = 0 when a < 0 . 
Note that when the sampling is uniform, that is, when tj = 7 , then Nnj{x) = Nn{x — j), 
where A„ is the cardinal B-spline of order n (361 Theorem 10.2.7] satisfying 


( 11 ) 



{x-k)l. 


To prove the theorem, note that the technique of using the Fourier transform of the fun¬ 
damental cardinal spline function does not work here, due to the non-uniform sampling. 
However, due to the condition of if/, we could control the difference between the non- 
uniform sampling case we discuss here and the uniform sampling case. Recall that the ^-th 
order spline interpolation for the non-uniform sampling is flnding the unique (36l Theorem 
10 . 2 . 2 (d)] ^-th order spline / such that 


N-n 

( 12 ) f(x) = Y, CjNnjix) 

j=-n 

where {c_„,,..., caz-^-i, CN-n} C M satisfles the conditions 


^ {U) — f{U ) 


j=-n 


(13) 
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for all / = 0, (361 (10.3.7)]. By the Schoenberg-Whitney Theorem (361 Theorem 
10.3.2], the interpolation is solved by finding cj by 

(14) c = A-^f, 

where A is a (A + 1) x (A + 1) matrix 

••• ^n,N—ni^o) 

••• ^n,N—n{H) 

A= ; : : : ’ 

J^n,—n{^N) ••• ^n,N—n{^N)_ 

C=[C-ri^C-ri^l,...,CN-n-UCN-nV andf= [/(^o),/(O /(W-n- 1 ),/(W-n)] ^ ^ 

Similarly, denote the ^-th order spline interpolation from the uniform sample 
to be /, which satisfies 


N-n 

(15) f { x )= ^ djNn { x - j ) 

j=-n 

where ..., C M satisfies the conditions 


N-n 

(16) ^ djNn{i-j) =f{i) 

j=-n 

for all / = 0,..., A. The dj is again solved by applying the Schoenberg-Whitney Theorem. 
By the assumption of if/, we know that \f{ti) — f{l) | < C|/ — T|e for all / G / by Q, where 
C' depends on ||/'||oo. Also, we have \Nnj{x) — Nn{x — j) \ < C"\j — t\£ = 0{^^) for all 
j G /, where C" is another constant depending on n. Hence, cj = dj + by the ^/£- 

perturbation of A. Also note that Nnj{x) and Nn{x) are all compactly supported. As a 
result, over [0,A], the ^-th order spline interpolation over the non-uniform sampling ^ 
satisfies 

fix) = Y! CjNn,jix) = djNnix-j) + OiVe) 

j=-n j=-n 

(17) =/(x) + 6>(Ve). 

Thus, by the above argument, we know that f{x) satisfies the reflection property ([^, and 
hence we obtain the proof. 

□ 


4. Reflection effect and time-frequency analysis 

In practice, the refiection effect discussed in Sectionj^is pronounced when we apply the 
TF analysis techniques to study the time-varying dynamics hidden inside the signal, and the 
effect is even worsened when we apply the sharpening technique in the TF analysis. In this 
section and the next, we demonstrate how the interpolation induced refiection effect is pro¬ 
nounced by two TF analysis techniques, the RM (571[T0l[TTl or the SST (l2l[T3l. Similar 
phenomena could be found in other TF analysis methods, ranging from the linear to qua¬ 
dratic methods 0. Note that the SST is a variation of the RM, and these techniques could 
be carried out to sharpen the TF representation determined by a chosen linear TF analy¬ 
sis, for example, the short time Fourier transform (STFT) or continuous wavelet transform 
(CWT). In a nutshell, RM and SST are nonlinear TF techniques aiming to alleviate the 
spreading effect in the TF representation determined by the linear TF analysis, which is 
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caused by the inevitable Heisenberg uncertainty principle. These techniques sharpen the 
TF representation by taking the phase information hidden inside the linear TF analysis into 
account. For details, we refer the reader with interest to (381 [T3 for the review material. 



N 

I 



30 40 50 

Time (sec) 


10 20 


30 40 50 

Time (sec) 


60 70 



Time (sec) 


Frequency (Hz) 


Figure 2. The time-frequency (TF) representation of fi{t) = 
cos(2;r2.5r) and different interpolations from the sampling scheme y/i. 
The instantaneous frequency (IF) of fi is (j)' = 2.5. Left: the SST-STFT 
result of fi (t). Middle left: the SST-STFT result of the cubic spline (CS) 
interpolation. The instantaneous Nyquist frequency (INF), is su¬ 

perimposed as a red dashed curve on it. Here we could see components 
with IF \i/[ — 0' marked as (1) and i/rj + (j)' marked as (2). Middle right: 
the reassigned STFT result of the CS interpolation, where we could only 
see an extra component with IF — 0' marked as (1). Note that the ar¬ 
tificial component with IF — 0' is the refiection of (j)' associated with 
the INF. Right: from top to bottom we show the 40-th second slice of the 
TF representations shown in left, middle left, middle and middle right. 
The unit of the y-axis is arbitrary (a.u.). 


To demonstrate this refiection effect, we start from a harmonic function fi (t) = cos(2;r2.5r) 

and sample it with the ISR i/rj {t) = 6 • Note that obviously the ISR is greater 
than the INR of /. Then, we demonstrate the refiection effect by applying different com¬ 
mon interpolation schemes. We then run SST-STFT and reassigned STFT on the interpo¬ 
lated signal, where the interpolated signal is sampled uniformly at 64 Hz. To avoid possible 
boundary effects, we sample the signal for 80 seconds. The results are shown in Figure 
In this study, the TF representation R e is displayed in the log scale. Precisely, we 
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plot R e where Rij = max{10“^,log(l + min{|7?/j|,^)})}, where q is the 99.8% 

quantile of all entries of |/?|. 

We next consider the same procedure on a non-harmonic function f 2 {t) = (0.7 ) cos{27r{7rt + 

0.2cos(r))), whose AMandlFarea(r) =0.74-^^-^ and0'(r) = ;r —0.2sin(r), respectively, 
and we take another ISR, v4(0 = 8 + 0.5cos(;rr/10). Note that Xj/^ is greater than INR of 
/ 2 . The result is shown in Figure Note that we could see a clear reflected component 
associated with the INF in all the above cases, and the behavior depends on the setup. 




Time (sec) 



Time (sec) 


Figure 3. The time-frequency representation of f 2 {t) = (0.7 + 
r^-^)cos(2;r(;rr 4 -0 .2cos(r))) and the interpolation from the sampling 
dataset associated with V^2(0' The instantaneous frequency (IF) of /2 
is = 71 — 0.2sin(r). Left: f 2 {t) is shown as a gray curve with the 
non-uniform samples superimposed as black circles. The cubic spline 
(CS) interpolation is shown as the black curve. The unit of the y-axis is 
arbitrary (a.u.). Middle left: the SST-STFT result of /2(0- Middle right: 
the SST-STFT result of the CS interpolation. Right: the reassigned STFT 
result of the CS interpolation. Note that in the middle, middle right and 
right subflgures, in addition to we could see components with IF 
1/4 — marked as (1), 2i/4 — 0' marked as (2) and 1/4 + 0' marked as 
(3). Note that the artiflcial component with IF 1/4 — 0' is the reflection 
of 0' associated with the INF. 


5. Real signal from Anesthesia 

General anesthesia is usually inevitable for a patient receiving major surgery. For 
the short term and long term well-being of the patient, the anesthetic agent dose should 
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be dynamically adjusted to achieve an adequate level of anesthesia. It has been shown 
that the oscillatory patterns in the R-to-R peak interval (RRI) time series of electrocar¬ 
diography and the respiration contain a lot of information regarding anesthesia dynamics 
It has been shown in EHinilloliiel that the adaptive harmonic 
model, the multi-taper reassigned STFT |[l4l|4Tl and the multi-taper SST-STFT ifTSlfT^ 
serve as a good framework toward the goal, especially when noise is inevitable. Precisely, 
the oscillatory behavior of the RRI during anesthesia could be modeled by the adaptive 
harmonic model, and how strong the component is related directly to the anesthetic depth. 
We would extract stable features by the multi-taper reassigned STFT and the multi-taper 
SST-STFT to better quantify the anesthetic depth. While its clinical value has been shown 
in different problems, its sampling theory issue is left unanswered to the best of our knowl¬ 
edge. We mention that the sampling issue for the power spectrum approach to study HRV 
based on the stationarity assumption has been widely studied; see for example 14^ |4^ . 
It could be argued that reassigned STFT and SST-STFT add complications on the top of 
a spectrum approach, but note the difference between the underlying models which are 
aiming to capture different phenomena. 

Here we demonstrate two examples regarding this direction which might generate po¬ 
tential artifact in the TF representation - IHR analysis and FDR analysis. The IHR and 
FDR are acquired from the recorded ECG signal in the following way. Denote the recorded 
lead II ECG signal as E{t) which is digitalized at the sampling rate lOOOHz. The R peak 
detection was automatically determined from E{t). The ECG signal is clean without sig¬ 
nificant noise contamination, and no ectopic beats nor electro-cauterization happen in the 
recorded signal. The collected RRI time series is denoted as ^ where 

ti G M is the time stamp of the i-th R peak. Then, we follow the common practice and 
approximate the IHR from ^ by the cubic spline interpolation (21, and denote the approx¬ 
imated IHR as fm. Next, fm is resampled to be equally spaced at 8 Hz for the multi-taper 
SST-STFT and multi-taper reassigned STFT analysis (431 . as it is commonly believed that 
most useful information inside IHR is below 0.5 Hz under the stationary assumption. In 
addition to the time stamps of the R peaks, we also have a non-uniform sampling dataset 
{ti^E{ti)}f^^, where E{ti) is the amplitude of the i-th R peak. The FDR signal, a surrogate 
of the respiratory signal denoted as R(r), is built up by applying the cubic spline interpo¬ 
lation on {ti,E{ti)}f^^. We mention that although the amplitude scales of R{t) and R{t) 
are different, they share the same oscillatory information inside the respiratory signal, in 
particular the IF. We refer readers with interest in FDR to (25l|23| for details. To confirm 
the existence of the reflection effect as an artifact in the EDR signal, when we record the 
ECG signal, we simultaneously record the airway fiow signal, which is denoted as R{t). To 
avoid any possible artifacts, the airway fiow signal is uniformly sampled at the sampling 
rate 25 Hz. As R peaks are viewed as the non-uniform sampling of the IHR and EDR, the 
ISR, denoted as could be estimated by the cubic spline interpolated function from 

As we have the airflow signal serving as the ground truth for the respiratory signal, 
we start from showing the result of the EDR signal. The results of multi-taper SST-STET 
and the multi-taper reassigned STET on R{t) based on the cubic spline interpolation and 
the airway fiow signal are shown in Eigure Here the multiples of the base oscillatory 
component with the IE at about 0.5Hz are associated with the notion called “wave-shape 
function”, for which we refer reader with interest to (441. 

The analysis results of the IHR signal based on multi-taper SST-STET and multi-taper 
reassigned STET are shown in Eigures[^and|^ In Eigurethe whole 30 minutes analysis 
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Figure 4. The results of the ECG-derived respiration (EDR) signal and 
the airflow signal recorded simultaneously. The signals here are the same 
as those shown in Eigure while we show here the whole 280 second 
signal. The instantaneous Nyquist frequency (INE) is superimposed as 
a red dashed curve. Left: the multi-taper synchrosqueezed STET (SST- 
STET) result of the EDR based on the cubic spline (CS) interpolation. 

Right: the multi-taper SST-STET of the airflow signal. It is clear that 
the base component with IE about 0.5 Hz in the airflow signal is well 
captured by the EDR signal, while in the EDR signal there is an artifi¬ 
cial reflected component associated with the INE. Also note the temporal 
reassignment effect in the multi-taper reassigned STET around 100-th 
second in in Eigurethe TE resolution of the multi-taper reassignment 
STET is sharper than that of the multi-taper SST-STET shown, so the 
reflection effect is more visible. 

result is shown, while in Eigurej^we illustrate a zoom-in results of a 6 minutes sub-interval. 
The TE representation provides several informations, in particular the time-varying fre¬ 
quency of the IHR signal. Note that there is a dominant curve near 0.5 Hz, denoted as vi, 
which is associated with the well-known phenomena of respiratory sinus arrhythmia (the 
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respiratory signal is not shown), which oscillates at frequency about 0.5 Hz. It is clear 
that we can see the reflected effect in the interpolated IHR fm. Indeed, the dominant curve 
around 0.7 Hz is the reflection of the component xi which is related to the INF i/A'(r)/2. 



800 1000 
Time (sec) 



Time (sec) 


Figure 5. The time-frequency representation of the IHR signal, 
over the whole 30 minutes period. The instantaneous Nyquist frequency 
is superimposed as a red dashed curve on it; the IHR signal is super¬ 
imposed as a blue curve. Top: the multi-taper SST-STFT result of the 
IHR based on the cubic spline (CS) interpolation; bottom: the multi¬ 
taper reassigned-STFT result of the IHR based on the CS interpolation. 
The signal and the TF representation in the red box are zoomed in and 
displayed in Figure 


To close this section, we demonstrate one more example with a different interpolation 
scheme, the shape-preserving piecewise cubic interpolation (PCHIP). The TF representa¬ 
tion of the FDR signal generated via the PCHIP interpolation scheme determined by the 
multitaper SST is shown in Figure |7] It is clear that in the PCHIP interpolation scheme, 
the reflection effect still exists. This demonstration shows that the reflection effect is not 
specialized to the spline interpolation. In Figure [7] we also show the TF representation 
of the FDR signal generated via the PCHIP interpolation scheme determined by STFT. 
The STFT is carried out with the Gaussian window. It is also clear that we could see the 
reflection effect. 


6. Possible solutions to the reflection effect 


A naive solution to the reflection effect is by the hard threshold. Precisely, given a TF 
representation R:Rx M+ ^ C, we could deflne a new TF representation R by setting the 
TF representation coefficients with frequency larger than INF to zero; that is. 


(18) 


6/, eN / when ^ < w'{t)/2 

^ \ 0 when I > V/'(0/2 ' 


This will directly remove the reflection artifact induced by the spline interpolation. Al¬ 
though such an adaptive Altering could mitigate the problem, due to the nonstationarity, 
this approach might also remove the possible information hidden inside the signal. Another 
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Figure 6. The time-frequency (TF) representation of the IHR signal, 
fm(0. over a 6 minutes period. This is a zoom-in illustration of Fig¬ 
ure for the sake of showing the relationship between the IHR signal 
and the time-frequency representation. The instantaneous Nyquist fre¬ 
quency is superimposed as a red dashed curve on it; the IHR signal is 
superimposed as a blue curve. Top: the multi-taper SST-STFT result of 
the IHR based on the cubic spline (CS) interpolation; bottom: the multi¬ 
taper reassigned-STFT result of the IHR based on the CS interpolation. 
Note that the temporal reassignment in the multi-taper reassigned-STFT 
sharpen the reflected component around the period around the 1350-th 
second. It is clear that the IHR oscillates faster after 1350-sec, which 
leads to a higher instantaneous frequency. 


naive solution is to take a higher order spline interpolation so that in the Fourier domain the 
spectrum of the fundamental cardinal spline is closer to the step function. Indeed, it is well 
known that when n ^ oo^ f\^^) converges to the sense (271 (16)]. See Figure 

[^for an example of the FDR signal generate(i by the 12-th order spline interpolation. It 
is clear that compared with the cubic spline interpolation result shown in Figures and 
the reflection effect is alleviated. While this approach seems to work, however, it is well 
known that the higher the order of spline interpolation is, the more severe the overfltting 
is. This fact might limit its application. Depending on the application, we could consider 
different penalty to determine the optimal order of spline interpolation. Yet another naive 
solution is to pre-process the signal by a low pass Alter to the interpolated signal to reduce 
















WHEN INTERPOLATION-INDUCED REELECTION ARTIEACT MEETS TIME-EREQUENCY ANALYSIS 


15 


2 


^ 1.5 

N 

E. 


3 

O" 

<D 

-0.B 


0 



50 100 150 200 250 

Time (sec) 



50 100 150 200 250 

Time (sec) 


Figure 7. The time-frequency representation of the ECG-derived res¬ 
piration (EDR) signal generated via the piecewise cubic interpolation 
(PCHIP) interpolation scheme. The underlying ECG signal for the EDR 
signal is the same as that shown in Eigurej^ The instantaneous Nyquist 
frequency (INE) is superimposed as a red dashed curve. Left: the multi¬ 
taper reassigned short time Eourier transform (STET) result of the EDR 
based on the PCHIP interpolation. Right: the STET of the EDR based 
on the PCHIP interpolation. It is clear that the reflection effect exists. 
Also note that the multi-taper reassigned STET is sharper than the STET. 


the reflection artifact. However, this approach only works when the signal is band-limited 
- in general, for example the adaptive harmonic model, the structured artifact might not 
be removed but perturbed by the low pass filter, which leads to more complicated artifacts. 
These findings suggest that we should consider to design a different interpolation scheme 
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Figure 8 . The time-frequency representation of the ECG-derived res¬ 
piratory (EDR) signal determined by the 12-th order spline interpolation. 
The ECG signal we use to construct the EDR signal here is the same 
as that used for the EDR signal shown in Eigure The instantaneous 
Nyquist frequency is superimposed as a red dashed curve on it. Left: 
the multi-taper reassigned-STET result of the EDR signal determined by 
the 8-th order spline interpolation; right: the multi-taper SST-STET re¬ 
sult of the EDR signal determined by the 8-th order spline interpolation. 
Clearly, compared with Eigures[^and|^ the reflection effect is reduced. 


to balance between the the reflection effect and the interpolation purpose. This opens a 
future research direction in the TE analysis. 
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7. Discussion 

While extracting oscillatory features from a given time series via TF analysis is at¬ 
tracting more interest in bio-medical field, in particular when people want to study the 
dynamics, the artifacts are enhanced, in particular the reflection effect we discuss in this 
paper. Thus, when apply the TF analysis to study the bio-medical signal, in addition to 
search for features inside the TF representation, it is important to pay attention not only 
to the sampling effect but also to the interpolation scheme, which might generate artifacts 
and introduce misleading interpretations. 

The evidence of the reflection effect is shown clearly in the simulated results in Section 

In addition, it is clear that a reflection component exists in the FDR signal, as is shown 
in Figurej^in Sectionj^ while the reflection component does not exist in the airflow signal. 
While the reflection component has a quite different IF behavior compared with the base 
component, if not being careful, these observations might lead to a misinterpretation that 
there is a possible hidden structure inside the respiratory signal. The same comment holds 
for the IHR signal. The reflective component shown in Figure might lead to a misinter¬ 
pretation that there are two oscillatory components inside the IHR; however, it is actually 
originated from the numerical interpretation. In order to get the correct information, few 
possible methods to alleviate the reflection effect is proposed in Section 

In addition to designing a different interpolation scheme raised in Section there are 
several open problems and their importances have been indicated in this paper. First, while 
there are several different interpolation schemes available for reconstructing the signal 
from the (non-)uniform sampling points, depending on the application, it is important to 
pay attention to the interplay between the TF analysis and the artifact generated by the in¬ 
terpolation scheme, so that the results are not misled by these artifacts. Second, it raises a 
question specific to the HRV or BRV - as the sampling scheme is limited by the physiolog¬ 
ical facts, it is not possible to sample the instantaneous heart rate as fast as possible. Thus, 
is it possible to evaluated the dynamical spectral information directly from the sampled 
time series, so that we could avoid the need of the interpolation? One possible direction is 
the following. For a given non-uniform sampled time series, the Lomb periodogram ll4^ 
or the spectrum of counts EH are commonly applied methods if we want to estimate the 
power spectrum directly from the non-uniform sampled time series. While they allow us 
to avoid the interpolation, however, the necessary phase information for the reassignment 
technique and SST is missed in these methods. Thus, it deserves a further study to explore 
the possibility to combine the benefits of the Lomb periodogram or the spectrum of counts 
and the reassignment technique. 


8. Conclusion 

We report a potential artifact shown in the interaction of the interpolation and the TF 
analysis. This artifact is theoretically studied under the spline interpolation scheme, and 
is referred to as the reflection effect. A solution as well as the future directions are pro¬ 
vided. While in this paper we demonstrate the evidence based on the bio-medical signals, 
including the IHR signal and the EDR signal, this theoretical phenomena might happen in 
other fields. In conclusion, when we apply the TF analysis techniques to study the time- 
varying dynamics hidden inside a time series, it is important to pay attention to avoid any 
mis-interpretation, in particular when the sampling rate is not significantly higher than the 
spectrum we are interested in. 
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